{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial4: Convolutional Layers - Spectral methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Why convolution in ML\n",
    "- Some theory on convolution\n",
    "- Convolution on graphs\n",
    "- Spectral-convolutional layers in PyTorch Geometric"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.11.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/antonio/anaconda3/envs/geometric_new/lib/python3.9/site-packages/tqdm/auto.py:22: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n",
      "  from .autonotebook import tqdm as notebook_tqdm\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import torch\n",
    "os.environ['TORCH'] = torch.__version__\n",
    "print(torch.__version__)\n",
    "\n",
    "!pip install -q torch-scatter -f https://data.pyg.org/whl/torch-${TORCH}.html\n",
    "!pip install -q torch-sparse -f https://data.pyg.org/whl/torch-${TORCH}.html\n",
    "!pip install -q git+https://github.com/pyg-team/pytorch_geometric.git"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams.update({'font.size': 16})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why convolution in ML"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Weight sharing\n",
    "- Detection of translational invariant and local features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](fig/fully_connected.png)\n",
    "[Source](https://missinglink.ai/guides/convolutional-neural-networks/fully-connected-layers-convolutional-neural-networks-complete-guide/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](fig/cnn.gif)\n",
    "[Source](https://commons.wikimedia.org/wiki/File:Convolutional_Neural_Network.gif)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](fig/Convolution.png)\n",
    "[Source](https://commons.wikimedia.org/wiki/File:Convolution.PNG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some theory on convolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](fig/Convolution_of_box_signal_with_itself2.gif)\n",
    "[Source](https://en.wikipedia.org/wiki/File:Convolution_of_box_signal_with_itself2.gif)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "c[n] = (v * w)[n] = \\sum_{m=0}^{N-1} v[m] \\cdot w[n-m]\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def conv(v, w):\n",
    "    c = np.zeros(v.shape)\n",
    "    for n in range(len(v)):\n",
    "        c[n] = 0\n",
    "        for m in range(len(v)):\n",
    "            c[n] += v[m] * w[n - m]  \n",
    "    return c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD9CAYAAACcJ53WAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAzEklEQVR4nO3deXRcV5nv/e8uVWmybMmyLMnWUI7jeHY8KZ4HSS45cQCnuUC4zSWwLjeE5gXe5vYLNDT07e7btxto4NIsFs0iDms1nYQh0EBMp+2kSpPneLbjOZ5KgzXYliVrlkq13z/KZRzHiqqOqurUqXo+a3nJLmnX+fm4/Ojoqb33UVprhBBCWIfN7ABCCCHCI4VbCCEsRgq3EEJYjBRuIYSwGCncQghhMfZoHyAvL0/PmDHD8Pje3l4mTJgQuUBJRs7f+Mj5Gx85f8YdOXLkhtZ66oM+F/XCPWPGDA4fPmx4fF1dHeXl5ZELlGTk/I2PnL/xkfNnnFLKO9rnpFUihBAWI4VbCCEsRgq3EEJYjKHCrZTaqZTSSqn/E+lAQggh3lvYhVsp9afA4ihkEUIIEYKwCrdSKgf4PvAXUUkjhBBiTOFOB/wn4LTW+hdKqZ9HI5AQIuB4+3EOtx2mrKCMJflLzI5j2ODgIB0dHXR3dzMyMmJ2HNOlpqaSl5dHdna24ecIuXArpdYBn0DaJEJE3fH243zq9U8x4h8hNSWVbZu3WbJ4K6VoaGhg8uTJzJgxA4fDgVLK7Fim0VrT399PU1MTaWlppKenG3qekAq3UsoB/AT4rtb6fAhf/xzwHEBBQQF1dXWGwgH09PSMa3yyk/M3Pmadv1dvvcqwfxiAoZEhXtn/Cp3ZnTHPMV42m4309HTS0tIYHBxkcHDQ7EhxYcKECZw8eZK+vj5D40O94v5LIAP4h1C+WGv9PPA8QFlZmR7PyilZeTU+cv7Gx6zzd/DgQbgd+L0jxcHTq5+25BX30aNHKSgoIDU11ewocSU9PZ2+vj5WrFhhaPyYhVspVQp8HXgWSFNKpd3z6bQ7b1h2a62leSVEhJztOItd2fFpH59f+nlLFm0ItEocDofZMeKO3W7H5/MZHh/KrJKZQDrwEnDrnl8AX7rz+0WGEwgh3uFm/02Oth/lEws+QZYji8udl82ONC7J3NMezXjPSSitkuNAxQMeryVQzH8KXBxXCiHEXTWNNfi1nycfepK2vjZqG2vx+X3YbVHfE05YxJivBK11J1B3/+N3vmN4tdbv+pwQwjiP10PJxBJmT55NVWkVr11+jcNth1k1bZXZ0USckL1KhIgjXYNdHGw5iMvpQinFmqI1ZNgz8Hg9ZkcTccRw4dZaK631NyIZRohkV9dYh0/7qCqtAiDDnsG6onVUN1Tj135zw4m4IVfcQsQRj9dD4YRCFuYtvPtYlbOKG/03ON5+3LxgIq5I4RYiTvQO97Lv2j5cpa53zDrYULyBVFsqbq/bxHTx5Yj3Fj+qvcgR762xvziCXnnlFZRSnDx58l2f27JlC0uWLIlJDnmbWog4satpF0P+IVxO1zsen+CYwJrpa6huqOYrj30lIabX/d0fTnPm2m1DY7sHhjnX2o1fg03B3MKJTEwPf674/OmT+JsPLAhrzNatW8nOzuall17in/7pn+4+3tbWhsfj4Vvf+lbYOYyQK24h4oTb62ZK+hSWTF3yrs+5nC5aels4ffN07IPFmdsDPvw68Hu/Dvw5VtLT0/nIRz7Cz3/+c/z+P77n8Itf/AKtNR/72MdikkOuuIWIA/2+fvY07+EDMz9Aii3lXZ8vLynHruy4ve539L+tKtwr3Xsd8d7iv71wgGGfH4fdxg/+61KWOydHMN17e+aZZ3jhhReoqanB5Qr8dPTiiy/icrmYNm1aTDLIFbcQcWBf8z76ff3vapMEZadls2LaCjxeD1rrGKeLL8udk3n52VX8xeY5vPzsqpgWbYD169czY8YMXnzxRQDOnj3L0aNHeeaZZ2KWQQq3EHHA3eAmOy2bssKyUb/G5XTR0N3AhVsXYpgsPi13TuZzFbNiXrQhsPjw4x//OL/97W/p6+vjxRdfJCsriw9+8IMxyyCFWwiTDY0MUd9YT0VJBQ7b6G+yVZZUYlM2PA2yGMdszzzzDD09Pfz2t7/l5Zdf5kMf+hCZmZkxO74UbiFMdqDlAD3DPVQ5q97z66ZkTGFZ/jJZRRkHZs+ezcqVK/nqV79KQ0NDTNskIIVbCNN5vB6yHFkh7UXicrq42HmRK11XYpBMvJdnnnmG5uZmioqKqKh40D580SOFWwgT+fw+ahtrA4tsUsa+2YCrNPDmpVx1m+9zn/scWmuampqw2WJbSqVwC2Giw22H6RzsHLNNElQwoYBHpz4qqyiTnBRuIUzk8XrIsGewtmhtyGOqSqs423GWpu6mKCYT8UwKtxAm8Ws/1Q3VrCtaR4Y9I+Rxwbne1Q3V0Yom4pwUbiFMcrz9ODf6b9ztW4eqeGIx83LnSbskiUnhFsIkbq8bh83BhuINYY91OV2cuH6Ctt62KCQT8U4KtxAm0FpT3VDNmulryErNCnu8tEuSmxRuIUxw+uZpWnpbRt2bZCwzs2fycPbDsooySUnhFsIEbq8bu7JTUWJ84YbL6eJI2xE6BjoimExYgRRuIWJMa43H6+GxwsfITss2/DxVzir82k9NQ00E0wkrkMItRIxduHWBhu4Gw22SoNmTZ1MysURWUSYhKdxCxJinwYNCUVlaOa7nUUrhcrp4s+VNuga7IpROWIEUbiFizOP1sKxgGXkZeeN+rqrSKnzaR31TfQSSCauQwi1EDF3pusLFzosh700yloV5CymcUCiLcZKMFG4hYijYj95Uuikiz6eUwlXqYl/zPnqHeyPynJbQeBB2fy/wMYYOHz6MUoo9e/bcfeyHP/whSim+8Y1v3H3s7bffRinFf/7nf0Ylh9wsWIgYcnvdPJr3KIUTCiP2nC6ni5fOvsTupt088dATEXveqNrxVWh9y9jYwdvQdgq0H5QNChZC2qTwn6dwEWz5VlhDli1bRk5ODjU1Naxbtw6AmpoaMjIyqKn54+yempoaUlJSWL9+ffi5QiBX3ELESFN3E2c7zo57Nsn9lkxdwpT0KcnTLhnoChRtCHwciN0bszabjQ0bNlBbWwuA3++nvr6ez372sxw6dIienh4AamtrKSsrY+LEiVHJIVfcQsRIcHl6pAt3ii2FTaWb+MPlPzDgGyDdnh7R54+KMK9036HxIPxsK4wMQUoqfOgFKFkRuWxjqKio4Ktf/SoDAwOcOXOGzs5OvvKVr/CTn/yE3bt3s2XLFurq6vjUpz4VtQxyxS1EjLi9bubmzqVkYknEn9vldNHv62fvtb0Rf+64U7ICPrkdKr8e+BjDog1QWVnJ4OAg+/bto7a2lsWLF1NQUMC6deuora3l9OnTtLW1RfV2ZnLFLUQMtPW2ceL6CT6/5PNRef6ywjKy07LxeD0Re+MzrpWsiHnBDlq0aBF5eXnU1NRw7NgxKisD8/ErKyt55ZVXKCkpITU1lbVrQ785RrjkiluIGAi2SSI1DfB+DpuDipIK6hvrGR4ZjsoxRIBSio0bN+J2u9m9e/c7CvexY8f43e9+x8qVK8nMzIxaBincQsSAp8HDzOyZzMyZGbVjVDmr6B7u5kDLgagdQwRUVlZy8OBB+vr67s4cWbZsGZMmTaK2tjbqd32Xwi1ElHUMdHCk7UjUWxirpq1igmOCbPUaA8HCXFZWxqRJgamIwRkn934+WqTHLUSU1TTU4Nf+qLVJglJTUtlQvIGahhr+etVfY7fJf+9omTdvHlrrdz3+6quvxuT4csUtRJR5vB6KsoqYmzs36seqclbROdjJkbYjUT+WMI8UbiGiqGuwizdb3qTKWYVSKurHWzt9Lekp6cmzGCdJSeEWIorqm+rxaV/EF92MJtORybqidXfbMyIxjVm4lVKPK6VqlFKtSqlBpVSTUuoVpdT8WAQUwsrcXjf5mfksylsUs2O6nC6u91/nxPUTMTumiK1QrrhzgSPA54HNwNeABcABpZQzitmEsLTe4V72Ne/DVerCpmL3w+3G4o04bA5plySwMV9NWutfaK2/rLX+jda6Xmv9IvBfgInAh6OeUAiL2t20myH/UMzaJEFZqVmsnr6aam/1A2c+COszehlw885HWaIlxCh+feHXZNgzSFEpMT+2q9TFtd5r/OOb/8jx9uMxP76IrpALt1IqRSmVqpR6BPgJ0Ar8MmrJhLCwgy0HOdh6kH5fP59xfybmxTM/Mx+AX53/FZ9+49NSvBNMODP03wSW3/n9RaBSa93+oC9USj0HPAdQUFBAXV2d4YA9PT3jGp/s5PyNj9Hz97MbP7v7+6GRIV7Z/wqd2Z2RCzaGN7reAECjTTl+0KRJk+ju7o75ca1gYGDA8P/NcAr3M8AkYCbwJcCtlFqntb56/xdqrZ8HngcoKyvT5eXlhsIB1NXVMZ7xyU7O3/gYPX+/dP8SeiFFpeCwOXh69dMsyV8S8XyjyWnPYcfrO/D5fdht9pgfP+jYsWNRu5mA1aWnp7N06VJDY0NulWitz2qt39Ra/wLYBGQBXzV0VCES2PDIMCevn2Rd0To+v/TzbNu8LeZFc0n+Ev65/J8BeN/M95lStEX0GNrMQGvdqZS6CMyKcB4hLO9AywG6h7v56JyPUl5SblqOjSUbWZa/jLduGLy3o4hbhmaVKKUKgLnApcjGEcL6PA0eMu2ZrJ6+2uwoVDmruNh5katdV82OIiIolJWTv1NK/bVS6imlVIVS6jNAPeADvhf1hEJYiM/vo6ahho3FG0lLSTM7zt055Im21evx9uO88NYLps2WOXHiBB/84AeZMmUKGRkZzJkzh29+85sxO34orZIDwNPA/wekAo1AHfDNB70xKUQyO9J2hM7BzpgvuhlN4YRCFuUtwu118+yiZ82Oc9e3D36bcx3nDI3tGerh/K3zaDQKxZzJc8hKzQr7eebmzuUvV/xl2OMOHjxIeXk5s2bN4vvf/z7FxcW8/fbbnDx5MuznMmrMwq21/jbw7RhkEcLy3F436SnprCtaZ3aUu1xOF98/8n2ae5opyioyO864dQ93owmsCNVouoe7DRVuo770pS8xZcoUDhw4cPf2ZMHbl8WK7LQuRIT4tZ+ahhrWFq0l0xG9+w2Gq6q0iu8f+T4er4dPLvik2XEADF3pBh1vP86n3/g0w/5hHDYH31r/rZjNmunr62Pv3r18+ctfjuo9JccihVuICDlx/QTX+6/HTZskqGRSCXMmz4mrwj0eS/KXsG3zNg63HaasoCymUx1v3bqF3++nuLg4Zsd8ECncQkSI2+vGbrOzsXij2VHexeV08aPjP6K9r/3ucngrW5K/xJS56ZMnT8Zms9Hc3BzzY99LbqQgRARoran2VrN62mompsbfSsHg/S6rG6pNTmJtmZmZrFu3jpdeeon+/n7TckjhFiICztw8w7Xea1G/IbBRD+c8zEPZD+HxJta0QDN897vf5ebNm6xevZoXX3yR2tpafvrTn/KFL3whZhmkcAsRAW6vmxSVQkVJhdlRRuUqdXG47TAdAx1mR7G0xx57jL1791JSUsIXvvAFnnzySb7zne/EtO8tPW4hxklrjafBQ1lhGTnpOWbHGVWVs4ptb22jtqGWD83+kNlxLG3p0qX84Q9/MO34csUtxDi93fk23tteqkrjs00SNDd3LkVZRbgb5JZmVieFW4hx8ng9KBSbnJvMjvKelFJUOat4s+VNbg/dNjuOGAcp3EKMk9vrZmn+UvIy8syOMiaX04XP76O+sd7sKGIcpHALMQ5Xu65ysfNi3C26Gc2ivEXkZ+bLHeAtTgq3EOMQ3HXPVWqNwm1TNlylLvZd20ffcF9Mjil3mn+38Z4TKdxCjIPb62bhlIVMy5pmdpSQuZwuBkcG2dW8K+rH8vv9pi5UiVf9/f04HA7D46VwC2FQc08zZ26esUybJGhZ/jJy03Njshinr6+PpqYmOjo6GB4eTvqrb601fX19NDc3k59vfOsBmccthEHBwhevqyVHk2JLobK0ktcuv8aAb4B0e3rUjjUyMkJpaSnXr1/n5s2b+Hy+qB3LKhwOBwUFBUyaNMnwc0jhFsKg6oZqZk+eTemkUrOjhK2qtIrfXPgN+6/tp6I0uqs909PTKSkpieoxko20SoQw4HrfdY63H7dcmyTosWmPMSl1UsLd0ixZSOEWwoDqhmo0Ou5XS47GYXNQXlJObWMtwyPDZscRYZLCLYQBHq+HGZNm8HDOw2ZHMazKWUX3UDcHWw+aHUWESQq3EGG6NXCLw22HqXJWoZQyO45hq6evJtOeKYtxLEgKtxBhqm2sZUSPWLa/HZSWksbG4o2Bv49/xOw4IgxSuIUIk9vrpiiriHm588yOMm4up4uOgQ6Oth81O4oIgxRuIcJwe+g2B1oO4Cp1WbpNErSuaB3pKenSLrEYKdxChKG+sR6f32f5NklQpiOTtUVrqfZW49d+s+OIEEnhFiIMHq+H/Ix8Hp36qNlRIsbldNHe387J6yfNjiJCJIVbiBD1Dfex99peNjk3YVOJ819nY/FG7Da73EjYQhLn1SdElO1u3s3gyKDl9iYZy8TUiayethpPgyfpN4GyCincQoTI4/WQm57LsvxlZkeJuCpnFc09zZztOGt2FBECKdxChGBwZJBdTbuoKKkgxZZidpyIqyipIEWlSLvEIqRwCxGCfc376PP1JVybJCgnPYeywjLcXre0SyxACrcQIfA0eJiYOpEVhSvMjhI1VaVVXL19lUudl8yOIsYghVuIMQyPDFPbWEtFSQWOFOO3m4p3m5ybUCjcDbIYJ95J4RZiDAdbD9I91G2ZGwIblZeRx9L8pdLntgAp3EKMwe11k2nPZE3RGrOjRJ3L6eLCrQs03G4wO4p4D1K4hXgPfu2ntrGWDcUbSEtJMztO1AV/qpC9S+KbFG4h3sOlwUt0DHQkzN4kY5mWNY2FUxZKuyTOjVm4lVIfVkr9u1LKq5TqV0qdV0p9Uyk1MRYBhTDT8b7jpKWksb5ovdlRYsbldHHq5ilaelrMjiJGEcoV95eAEeCvgCeAHwOfBdxKJdCGDULcx6/9nOg7wdrpa8l0ZJodJ2aCc9XlRsLxyx7C13xAa339nj/XK6U6gJ8B5UBNNIIJYbaT10/SNdKVNG2SoNJJpcyePBuP18Mz858xO454gDGvmO8r2kGH7nwsimwcIeLHz8/+HIViSsYUs6PEnMvp4mj7UX5w9Accbz9udhxxH6Otjo13PsqONCIhHW8/zs6rO9Fo/rzmz5OueJVklQDw07d+yqff+HTS/f3jXSitkndQShUB/xvwaK0Pj/I1zwHPARQUFFBXV2c4YE9Pz7jGJzs5f8b8+uav0QT27BgaGeKV/a/Qmd1pbqgY2tu5FwCNHtffX15/0RFW4VZKZQGvAj7gv4/2dVrr54HnAcrKynR5ebnhgHV1dYxnfLKT82eMe7cbekChSE1J5enVT7Mkf4nZsWImpz2HHTt3MKJHcKQ4DP/95fUXHSG3SpRS6cB2YCbwuNa6KWqphDCR1pqTN04yP3c+7895P9s2b0uqog2wJH8Jf7vmbwF4Zv4zSff3j3chFW6llAP4d2AF8KTW+q2ophLCRJc6L3H19lU++MgH2Zy9OWmL1lMPP0VRVhHnO86bHUXcJ5QFODbgZWAT8JTW+kDUUwlhIneDG4ViU+kms6OYSimFq9TF/pb9dA91mx1H3COUK+4fAR8Bvgv0KqVW3fOrOLrxhIg9j9fD4qmLmZo51ewopnM5Xfj8Puqb6s2OIu4RSuHecufj14H99/16Nkq5hDBFw+0GLty6kHSLbkbz6NRHyc/Il71L4syYs0q01jNikEOIuBDcFU8Kd4BN2agsreT3F39P33BfUi39j2ey14gQ9/B4PcyfMp+iLFkUHFTlrGJgZIA9zXvMjiLukMItxB0tPS2cunkqYW8IbNSygmVMTpss7ZI4IoVbiDuCu+El+i3KwmW32aksraS+qZ7BkUGz4wikcAtxl8frYVbOLGZkzzA7StxxOV30+frYf22/2VEEUriFAOBG/w2OtR+TNskoVhauZKJjotzSLE5I4RYCqPZWo9Eym2QUjhQH5SXl1DXWMewfNjtO0pPCLQSB1ZLOSU4eyXnE7Chxy+V0cXvoNodaDo39xSKqpHCLpNc50Mnh1sO4Sl0opcyOE7fWTF9Dhj0Dd4O0S8wmhVskvdrGWkb0iPS3x5BuT2dD8QZqGmoY8Y+YHSepSeEWSc/tdTN9wnTmT5lvdpS453K66Bjo4Gj7UbOjJDUp3CKpdQ91c6DlAJucm6RNEoINRRtIS0mTxTgmk8Itktqupl0M+4elTRKiTEcma6avwdPgwa/9ZsdJWlK4RVLzeD1MzZjK4qmLzY5iGVXOKtr72nnrhtxPxSxSuEXS6hvuY0/zHipLK7Ep+a8Qqo0lG7Hb7NIuMZG8WkXS2nttLwMjA9ImCdOk1EmsnLYSt9eN1trsOElJCrdIWm6vm5y0HJYXLDc7iuVUlVbR3NPMuY5zZkdJSlK4RVIaGhliV9MuKksrsdvGvJ+IuE9FaQU2ZZO9S0wihVskpf3X9tM73CtbuBqUm55LWUHZ3a1wRWxJ4RZJye11M9ExkVXTVpkdxbJcThdXuq5wqfOS2VGSjhRukXSG/cPUNtaysWQjjhSH2XEsa1PpJgBpl5hACrdIOodaD3F76LZs4TpO+Zn5LJm6RKYFmkAKt0g6Hq+HDHsGa6evNTuK5bmcLs7fOk/j7UazoyQVKdwiqYz4R6huqGZ90XrS7elmx7G84E8tstVrbEnhFknlWPsxOgY6ZNFNhBRlFTF/ynxpl8SYFG6RVDwNHlJtqawvXm92lIRR5azirRtv0drbanaUpCGFWyQNv/bj8XpYU7SGCY4JZsdJGMG58HLVHTtSuEXSOHXjFG19bdImibAZ2TOYlTNLpgXGkBRukTQ8Xg92ZWdj8UazoyScKmcVx9qPcaP/htlRkoIUbpEUtNa4vW5WTltJdlq22XESjsvpQqOpaagxO0pSkMItksL5W+dp6mmSRTdR8kjOIzgnOaVdEiNSuEVScHvd2JSNytJKs6MkJKUUrlIXh1oP0TnQaXachCeFWyQFj9fD8oLl5Kbnmh0lYVU5qxjRI9Q21podJeFJ4RYJ73LnZS53XZYtXKNs/pT5TJ8wXbZ6jQEp3CLhBfuuwd3sRHQopdjk3MT+a/vpGeoxO05Ck8ItEp6nwcPiqYspmFBgdpSEV+WsYtg/TH1TvdlREpoUbpHQGrsbOddxThbdxMjiqYuZmjFVVlFGWUiFWylVrJT6oVJqv1KqTymllVIzopxNiHELFhBpk8RGcObOnuY99A33mR0nYYV6xT0LeBq4BeyOXhwhIsvj9TAvdx7FE4vNjpI0qpxVDIwMsO/aPrOjJKxQC/curXWB1vpJ4NfRDCREpLT2tnLyxklpk8TY8oLl5KTlyGKcKAqpcGut/dEOEhWNB2H39wIfk3G8xR3x3uJHtRc54r1laPx39m8DYGhgsuHj/8elIcPHH29+q7Lb7FSWVlLbUMuOzh0cbz9udqSEYzc7QNQ0HoR/fT+MDEOKHTb/A+TPC318+1l44+sw4jNx/DfA74OUVPjkdihZEfp4izvivcXHth1geMSP3Wbj6++bx+yCiSGPr7t6iNebAj8c/vj0P3K7O4vyGY+FPP5CWzf/8NpZhkf8bL98IOzjB8f7/H5S7TZefnYVy53GvoFY0czsmfSP9LOjawc1b9SwbfM2luQvMTtWwlBa6/AGKPUssA14SGt9dZSveQ54DqCgoGD5L3/5S8MBe3p6yMrKCnvcI+d/TFHLTsPHjSd+bFx96L/R4Pxw2GONnj+z/fr8IK9d8Rken5r3Oql5tSgFWiuGrm9m6GZFBBOGzgb8l0ccvP/hVFOOb4adnTt5res1AGzYeF/O+9icvdnkVNZSUVFxRGtd9qDPReWKW2v9PPA8QFlZmS4vLzf8XHV1dRgaf/3foAVQNrCN44rZ7zNv/I6vABqbPY2ZlZ9gpoErbsPnz2R7es/AlSvYFIauuLedPsShW4Gijbbz35e5DF9xO1LCP/6Ftm7+7g+n8WtIddj4U9djSXXFndOew46dO/BrP44UB0+vflquuCMoMVslWkPTISheCXMehxnrw28zPLQepi+Bq7vNG998BE7+Cj76UlK1SQBON9+maHIGH1tRyqqZU8Iuev96uYXJvXk8kvE4m2eu5aOPhnerstUPT2FhUTa/8BwyVHRXPzyF1q5+flx/mb9/amFSFW2AJflL+MLSL/CDoz/gi8u+KEU7whKzcF87Bl0NUP6XsPTjxp+nZMX4CuZ4x6/6Mzj5S+hJrnv53ewZ5M0rN/lcxSw+VzEr7PFdg10cbDnIJxZ8gv+5/H8azrHcOZnuh1MNF91n18/kJ7su472ZnPOZPz7v4/z42I+50nXF7CgJJzFXTp7dDioF5jxpdpLxmbYEckrhzHazk8SU+0wbfg1PLCw0NL6usQ6f9pk+DXBKVhorH5rCjlMtpuYwS7o9nQUZC6huqGbEP2J2nIQScuFWSn1YKfVhYPmdh7bceSy+7gOldaDQPbQeMi2+hadSMG8rXK6Fgdtmp4mZHadaKc3NZP60SYbGe7wepk2YxoIpCyKcLHxbFhVy6Xovb7d1mx3FFIszF3Nz4CbHrx83O0pCCeeK+9d3fv3ZnT//y50//12kQ41L+xnouBQoeIlg3lYYGYILr5udJCa6+ofZd+kGWxYWopQKe3zvcC/7ru1jU+kmQ+Mj7fEFgZ8adpxKrnZX0IKMBaTaUmXvkggLuXBrrdUov8qjmC98Z7YDCua+3+wkkVH8GEycBmdfNTtJTFSfbWN4RBtuk+xq2sWQf8j0NklQwaR0ljsnJ23hTrels6ZoDZ4GD+FOPRajS7we99ntULoaJibIFp42W+Cb0NseGOo1O03U7TjVyrTsdBYX5xga7/a6ycvIi6tZDFsWFnK25Tbem4n/7/cgVc4qWntbOXXjlNlREkZiFe4bFwOtkvkJ0iYJmr8VfP1wMbF/3Owd9LHrwnUeX1CIzRZ+m6Pf18+e5j1sKt2ETcXPSzvZ2yUbizdiV3bcDbJ3SaTEz6s7EoLthHkfMDdHpJWugcwpCT+7pPZ8O4M+P1sMtkn2Ne+j39cfd3dyL8nNZFFRdtIW7uy0bFZOW4nHK+2SSEmswn1mOxQth+wE28IzxQ5z3xd4g9I3aHaaqNlxqpW8rFTKZhibDeRucJOTlkNZwQNXCZvqiYWFnGjs5Fpnv9lRTOFyumjsbuTCrQtmR0kIiVO4b3mh5XjizCa537ynYKgbLiXmHbQHhkeoPdfO5gWFpBhokwyNDFHfWE9FSQV2W/ytKwv+FLEzSa+6K0srsSmbbPUaIYlTuM/+IfAx0frbQQ9tgLTswJuvCWjXhev0DY0YbpMcaDlAz3BP3LVJgmZOzWJOwcSkLdy56bksL1gu0wIjJIEK93YoWAS5M81OEh32VJizBc69FtiqNsHsPNVKdoaDVTOnGBrv8XrIcmSxatqqCCeLnCcWFnLI20F794DZUUzhKnVxqesSl7sumx3F8hKjcN9ugcY3E/dqO2j+VhjoDGxclUCGfH7cZ9uoml+AIyX8l6TP76O2sZaNJRtJTYnfrVO3LCpEa3jjdJvZUUwRvO+nXHWPX2IU7nP/EfiYqP3toIcrwTEh4WaX7Lt0g+4Bn+E2yeG2w3QOdlJVGh+LbkYzp2AiD+VNSNp2ScGEAhZPXSyFOwISo3CfeRXyZkP+XLOTRJcjA2ZvDnyjSqBNe3aeaiUrzc66R/IMjfd4PWTYM1hTtCbCySJLKcUTCwvZf/kmt3qHzI5jiipnFWc7ztLY3Wh2FEuzfuHuvQHevYl/tR00byv0XoeGA2YniQjfiJ83zrRROTefNHtK2OP92k91QzXritaRYc+IQsLI2rKwkBG/xn02udsl1d5qk5NYm/UL97nXQPsTv78d9MhmsKcnzOySg1c76OgdMtwmOd5+nBv9N+Jmb5KxLCrKpignI2nbJcUTi5mXO09WUY6T9Qv32e2Q44TCR81OEhtpWfDwpsD0R7/f7DTjtvNUK+kOGxvnTDU03u11k2pLZUPxhggni45gu2TP2zfoHki82UGhqHJWcfL6SVp7k/ObVyRYu3D3d8Ll+sDVdhxs4Rkz87fC7Wa4dtTsJOPi92t2nmqlfHY+manhL5rRWlPdUM2a6WuY4JgQhYTRsWVhIUMjfmrOtZsdxRTBufbVDdIuMcrahfvCTvAPB1YVJpPZTwRuQHzG2lu9Hmu8RXv3IFsWGWuTnL55mpbelrhddDOaZaWTmToxjR1vJecV50PZDzErZ5bMLhkHaxfuM9th4vTA/iTJJCMHHtoYaBNZeNOeHW+14khRVMzNNzTe7XVjV3bKS8ojGyzKbDbF4wsKqLvQTt+Qz+w4pnA5XRxtP8rN/ptmR7Ek6xbuwR64VB3YCdBm3b+GYfO3wq2r0PqW2UkM0Vqz41Qr62blMSndYWi8x+thxbQVZKdlRyFhdG1ZOI2BYT+7Llw3O4opXKUu/NpPTWON2VEsyboV7+03wDeQPLNJ7jf3/aBslp1dcqr5Ns2d/WxZOM3Q+Au3LtDQ3WC5NknQyodymZzpSNqtXmdPnk3pxFJplxhk3cJ9djtMmBq4200ympAHzrWWXUW541QLKTZF1XxjdyryNHiwKRuVJZURThYb9hQbVfMLqDnbzqAvcRZThUophcvp4mDLQboGu8yOYznWLNzD/XDhjcAe1bbwF20kjHlb4cZ5uH7e7CRh0Towm2TVzFwmTzC2t4jH62FZ/jKmZBjblCoebFk4je5BH3sv3jA7iimqnFX4tI+6xjqzo1iONQv3pRoY7k2e1ZKjmXfnhsgWu+q+0NbD5Ru9PGGwTXKl6woXOy9atk0StGbWFCam2ZN2dsmCKQsonFAo7RIDrFm4z2yH9JzAHtXJbNJ0KF5huTvA7zjVglLw+AKDbZI7/9GDy6etKs2ewqZ5+bjPtjE8Yv3FVOFSSuEqdbHv2j56h5PzRspGWa9w+4bg/A6Y8ySkhD8bIeHM3xqYWdJxxewkIdt5qpUy52TyJ6YbGu/2unk071EKJxib/x1Pnlg4jc6+Yd683GF2FFO4nC6G/EPsatpldhRLsV7hvrILBruSdzbJ/YI3RrbI7JIrN3o519ptuE3S1N3E2Y6zlm+TBG2cPZUMRwo7TrWYHcUUS6YuYUr6FLmlWZisV7jPvgqpWTCzwuwk8WHyDJi22DJ97mCBesLgplLBZdKJUrgzUlOomDuV10+3MeK37mIqo1JsKWwq3cSe5j30+5LzRspGWKtwj/gCuwHOfhwcxn7MTkjztkLzYehqNjvJmHaeamVxcWCHPCPcXjdzc+dSMrEkwsnM88TCadzoGeSI95bZUUzhcrro9/Wzr3mf2VEsw1qFu2Ef9N2U2ST3m39nr5bgDZPjVNOtPk42dRluk7T1tnHi+glcpYlxtR1UOTefVLstadslZYVlZKdly1avYbBW4T6zHewZ8Ig19l6OmbxHYOq8uO9zB/egNrr3drBNYpW9t0OVlWZnwyN5vH6qFW3hvWeMctgcVJRUUN9Yz9BIct4ZKFzWKdx+f+CKctYmSLXOFp4xM38rePdBT/xuFbrzVCtzCycyI8/Yv5+nwcPM7JnMzJkZ4WTme2LhNK51DXCiKTlXEVY5q+gZ7uFAS2Lc2SnarFO4mw5BT+sf2wLineZtBfQfb5wcZ9pvD3Ck4ZbhvUk6Bjo40nYkYd6UvF/VvALsNpW07ZJV01aR5ciSxTghsk7hPrsdbI7AG5Pi3QoWQO7MuJ1d8vrpVrTG8N7bNQ01+LU/4dokQdmZDlY/PIWdSdouSU0J3MWotrEWnz85t7oNhzUKt9aBgvRwBaRbbwvPmFAqcNV9dTf0xd9ijh2nWpk5dQKP5GcZGu/xeijOKmbO5DkRThY/tiychvdmH2dbus2OYooqZxWdg50cbjtsdpS4Z43C3XIcuhpkNslY5m8Fvy+wsjSOdPQO8eaVDrYsLEQZuMVc12AXb7a8SZWzytB4q9i8oACbgp1J2i5ZW7SWDHuGtEtCYI3CfWY7qJTAboBidNOXQXZJ3M0ucZ9pZcSvDfe365vq8Wlfwva3g/Ky0nhsRm7S7tGdYc9gXdE6qhuq8evk27slHPFfuLUOFKIZ6yAz1+w08U2pwBL4SzUwcNvsNHftONVK8eQMFkyfZGi82+umILOAhXkLI5ws/mxZWMjb7T1cbO8xO4opXKUubvTf4Hj7cbOjxLW4L9wTehvg5kXZmyRU87bCyFDgDkFxoKt/mL0Xbxhuk/QO97KveR8upwubivuX67gFFycla7tkQ/EGHDaH7F0yhrj/n5B3Yz+gYO4HzI5iDSUrIasgbu4AX3OujeERbXi15O6m3Qz5hxJuteRoCrPTWVqak7TtkqzULNZMX0N1Q3VSzq4JVUiFWylVopT6jVKqSyl1Wyn1W6VUabTDAUy9vg9KV8FEY3s3Jx2bLXA/yoseGOozOw073mqlYFIaS0tyDI13e93kpueyNH9pZIPFsS0LCzl97TYNN83/9zODy+mipbeF0zdPmx0lbo1ZuJVSmUANMBf4JPAM8AhQq5SK7hLGm5fI6vXKbJJwzd8Kw32B4m2i3kEf9Reu88SCQmy28NskA74BdjfvZlPpJlKS6BZ1wTdxd55OznZJRUkFdmWXdsl7COWK+9PATOBPtNa/11q/CmwFnMBnohnu7uyIedImCYtzHWTkmj67pO78dQZ9fsNtkr3X9tLv60/42ST3K8nNZMH0SUnbLslOy+axwsfweD3SLhlFKIV7K3BAa30x+IDW+gqwF4jq+vM3j/0bP5xcwKsnXjc0/lcnd/M/fv9tfnVyd3KNT7HD3Cc5evF1XjvzHV6tf8HQ8V+tf4H/9bOPGB6/Y89PKc//F9oaf2do/K+PbyNNpZB645Kh8TQehN3fC3w0OL7U+5txjTd6/C0LC1GNB6l74aucO2TsJ6dzhzzs/9lfmTq+79ivDI13OV00dDfwxX97yrTXr9nj34sa6zuaUqoVeFVr/Zn7Hv8X4CNa66nvNb6srEwfPhz+Sqj/+5uv8K89/4kGFOAYmYSf0DszIwzhT/njCkLbSC4phH5HcauPz6CTnpT+u+cv36dJ1aG3K4aUpt2uIjZ+2ghkhDErZABNsy3w2kzXmm09NpbYMkMez3A/dDZAMEFOKTjC2AP8zniNRo1jvNHj9/Z2k9F7DXUnQauaypAtLeTxqf5BCvV1y44/ljrE/yoMvN7j4fU7nvFpGv565hd5auOzIY8HUEod0VqXPehz9hDG5wIP2uG9A5g8ygGfA54DKCgooK6uLrSk97jSfhCdCSiF1posv48h8kMe30cbpASmNmsNyu8gPYnGO+hCp3D3/KVpRZ4v9MJxwx4o+pEan+qHIlvoN79o0v2BnweVYhjYa4PphD6PP3O4iQloFKDR9A5r+hzWGT8y1EsG+s6/v6afVNpt00Men+9vBguPv5DWHDh3cfL6Hc94H5o9Z35Ptp4V8vixhFK4IXDZcL9Rv/1orZ8HnofAFXd5eXnYwbrUx9l/+Z/xobFr+IvZz4X1HetXJ3fz90f+HI0PtJ2/Wvk3fPTR9Ukz/tX6F/j7e87fc7PD+44f6fHPhjn++Kmf8+lD/8gwGoeGtRu/Rv7Cj4U8nsaD8LPAnHaVkkrWf/0pWSUrwh7v9w1is6cZHm/0+OcOecj5jz/FoX0MY2f4fT9g9WOh9/rPHfIwaOHx7fUvkBpHr9/xjl83/08o31ge8vixhNIqaQN+H+tWCQT+8nvO/J518/8k7B8zIFD83ri8j80z14RV9BJl/HjP36v1L3Dk6ussn/G4KeOPn/o5hy+/TtnMx1kSTtEOajwY2HRrxnoIp+jeM/5yzb8xs/IThseP5/jnDnm4daaGyfMrmRtG0Yun8Zf3/Dsz133I0HizX39mj3+vVkkohbsGSNVar7vv8bo74ze+1/jxFG6Auro6jFyxiwA5f+Mj52985PwZ916FO5R3i7YDq5RSd287opSaAay98zkhhBAxFErh3gZcBV5VSj2llNoKvAo0Aj+JYjYhhBAPMGbh1lr3ApXABeBF4GXgClCptU7OLcyEEMJEIc0q0Vo3AB+KchYhhBAhiPvdAYUQQryTFG4hhLCYMacDjvsASl0HvON4ijzgRoTiJCM5f+Mj52985PwZ5xxtnUzUC/d4KaUOjzaXUYxNzt/4yPkbHzl/0SGtEiGEsBgp3EIIYTFWKNzPmx3A4uT8jY+cv/GR8xcFcd/jFkII8U5WuOIWQghxDyncQghhMXFXuJVSJUqp3yilupRSt5VSv1VKlZqdyyqUUuVKKf2AX51mZ4s3SqlipdQPlVL7lVJ9d87TjAd83WSl1AtKqRtKqV6llEcptciEyHEllPOnlJoxyutRK6VyzElufaHeAScmlFKZQA0wCHySwJ13/g9Qq5R69M6GVyI0/y9w6J4/+8wKEsdmAU8DR4DdwOb7v0AppQhsX/wQ8AUCt/H7GoHX5BKtdVPs4sadMc/fPb7Ju7eB7o5SroQXV4Ub+DQwE5gTvKu8Uuok8DbwGeD/mpjNas5qrQ+YHSLO7dJaFwAopZ7lwYVnK7COwG6YtXe+dj+BHTK/QuAbZLIK5fwFXZbXY+TEW6tkK3AgWLQBtNZXgL3AU6alEglJa+0P4cu2AteCRfvOuC7gDyT5azLE8yeiIN4K9wLg1AMePw3Mj3EWq3tZKTWilLqplPq5vE9g2Hu9JkuVUlkxzmNV31RK+e68d7Vd3iMYn3hrleQS6CHerwOYHOMsVtUFfA+oB24DS4G/AvYrpZZqrdvNDGdBuQTuAHW/jjsfJwNyQ5HRDRK4U9YbwHVgLoHX4z6l1Aqt9Vkzw1lVvBVuCLwheT8V8xQWpbU+Bhy756F6pdQu4CCBfuw3TAlmXQp5TRqmtW4B/uyeh3YrpXYS+Inl68DHTQlmcfHWKrlF4ArnfpN58JW4CIHW+iiBW889ZnYWC+pg9NckyOsybFrrRmAP8no0LN4K92kCPcX7zQfOxDhLohntylG8t/d6TTbIfVcNk9fjOMRb4d4OrFJKzQw+cGdC/1rePQdUhEgpVQbMBt40O4sFbQeKlFIbgw8opSYBH0Bek4bceaN8LfJ6NCyuNplSSk0ATgD9BHqxGvh7YCLwqFzdjE0p9TKBOcZHgU4Cb05+DegDlmmt5W4k91BKffjObzcR6MX+PwTeRLuuta5XStkI/FhfAnyZPy7AeRRYfOfH/qQVwvn7HoELxP13Hp9D4PxlAyu11udjn9r64qpww93vxt8Hqgj8OFUNfFFrfdXMXFahlPoa8KeAE8gEWoEdwN/ceaNI3EMpNdp/gHqtdfmdr8kFvgv8CZBOoAj9hdb6RCwyxrOxzp9S6lPAZwmsspxI4DZmNcDfSdE2Lu4KtxBCiPcWbz1uIYQQY5DCLYQQFiOFWwghLEYKtxBCWIwUbiGEsBgp3EIIYTFSuIUQwmKkcAshhMX8//q1DRGixnHlAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 20\n",
    "v = np.zeros(N)\n",
    "v[8:12] = 1\n",
    "w = np.zeros(N)\n",
    "w[1:5] = 1\n",
    "c = conv(v, w)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.gca()\n",
    "ax.plot(v, '.-')\n",
    "ax.plot(w, '.-')\n",
    "ax.plot(c, '.-')\n",
    "ax.legend(['v', 'w', 'c'])\n",
    "ax.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fourier transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transformation $\\mathcal F: \\mathbb{R}^N \\to \\mathbb{R}^N$ with\n",
    "\n",
    "\\begin{align*}\n",
    "\\mathcal F^{-1}(\\mathcal F (v)) &= v\\\\\n",
    "\\mathcal F(v * w) &= \\mathcal F(v) \\cdot \\mathcal F(w).\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This implies\n",
    "\\begin{align*}\n",
    "v * w &= \\mathcal F^{-1}(\\mathcal F (v * w))\\\\\n",
    "&= \\mathcal F^{-1}(\\mathcal F(v) \\cdot \\mathcal F(w))\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4.75972857, 3.84267706, 4.75547361, 4.69473343, 4.41980675,\n",
       "       4.38291796, 3.98819354, 4.36896188, 4.0910097 , 4.60256767,\n",
       "       4.9724941 , 5.12207561, 5.60658149, 4.99165311, 4.75559415,\n",
       "       4.99484192, 4.8482739 , 4.82087932, 4.5831123 , 4.63790021])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v, w = np.random.rand(N), np.random.rand(N)\n",
    "conv(v, w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4.75972857, 3.84267706, 4.75547361, 4.69473343, 4.41980675,\n",
       "       4.38291796, 3.98819354, 4.36896188, 4.0910097 , 4.60256767,\n",
       "       4.9724941 , 5.12207561, 5.60658149, 4.99165311, 4.75559415,\n",
       "       4.99484192, 4.8482739 , 4.82087932, 4.5831123 , 4.63790021])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.fft import fft, ifft # Fast Fourier Transform / Inverse FFT\n",
    "np.abs(ifft(fft(v) * fft(w)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition of the Fourier transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Fourier transform can be computed as\n",
    "\n",
    "\\begin{align*}\n",
    "\\mathcal F(v) = U\\cdot v, \\;\\;\\mathcal F^{-1}(v) = \\frac{1}{N}\\ U^H \\cdot v\n",
    "\\end{align*}\n",
    "\n",
    "where the $N\\times N$ matrix $U$ is defined as\n",
    "\\begin{align*}\n",
    "\\\\\n",
    "U = \n",
    "\\begin{bmatrix}\n",
    "u_0(0) & u_1(0) & \\dots & u_{N-1}(0)\\\\\n",
    "u_0(1) & u_1(1) & \\dots & u_{N-1}(1)\\\\\n",
    "\\vdots & \\vdots& & \\vdots\\\\\n",
    "u_0(N-1) & u_1(N-1) & \\dots & u_{N-1}(N-1)\\\\\n",
    "\\end{bmatrix} \n",
    "\\end{align*}\n",
    "\n",
    "and $u_0, \\dots, u_{N-1}$ are functions defined as\n",
    "\n",
    "\\begin{align*}\n",
    "u_n(x)&:= \\cos\\left(2 \\pi \\frac{n}{N} x\\right) - i \\sin\\left(2 \\pi \\frac{n}{N} x\\right).\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def matrix_U(N):\n",
    "    u = lambda n, N: np.cos(2 * np.pi / N * n * np.arange(N)) - 1j * np.sin(2 * np.pi / N * n * np.arange(N))\n",
    "    U = np.empty((N, 0))\n",
    "    for n in range(N):\n",
    "        U = np.c_[U, u(n, N)]\n",
    "    return U\n",
    "\n",
    "\n",
    "def fourier_transform(v):\n",
    "    N = v.shape[0]\n",
    "    U = matrix_U(N)\n",
    "    return U @ v\n",
    "\n",
    "\n",
    "def inverse_fourier_transform(v):\n",
    "    N = v.shape[0]\n",
    "    U = matrix_U(N)\n",
    "    return (U.conj().transpose() @ v) / N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.00000000e+00-0.00000000e+00j,  0.00000000e+00-2.22044605e-16j,\n",
       "        7.77156117e-16+0.00000000e+00j,  2.22044605e-16+1.05471187e-15j,\n",
       "        2.22044605e-16+6.10622664e-16j,  1.99840144e-15-6.66133815e-16j,\n",
       "       -1.11022302e-15-4.44089210e-16j, -2.83106871e-15-1.11022302e-16j,\n",
       "       -4.44089210e-16+1.60982339e-15j, -1.33226763e-15+4.44089210e-16j,\n",
       "        8.88178420e-16+3.88387152e-16j,  2.88657986e-15-3.55271368e-15j,\n",
       "        4.44089210e-16+4.99600361e-16j, -4.27435864e-15-1.66533454e-15j,\n",
       "       -2.66453526e-15+3.44169138e-15j,  6.21724894e-15-1.11022302e-15j,\n",
       "       -6.66133815e-16+1.60982339e-15j,  2.22044605e-16-1.44328993e-15j,\n",
       "        1.44328993e-15+4.32986980e-15j, -9.99200722e-16+2.22044605e-16j])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fft(v) - fourier_transform(v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 5.55111512e-17-0.00000000e+00j,  0.00000000e+00-4.16333634e-17j,\n",
       "        4.51028104e-17-1.38777878e-17j,  1.38777878e-17-5.55111512e-17j,\n",
       "       -1.73472348e-18-1.73472348e-17j,  8.32667268e-17+1.38777878e-17j,\n",
       "       -5.55111512e-17+2.42861287e-17j, -1.38777878e-16+0.00000000e+00j,\n",
       "       -5.55111512e-17-8.67361738e-17j, -5.55111512e-17-2.77555756e-17j,\n",
       "        6.93889390e-18-1.94193576e-17j,  1.45716772e-16+1.87350135e-16j,\n",
       "        0.00000000e+00-3.12250226e-17j, -2.27248775e-16+9.71445147e-17j,\n",
       "       -1.24900090e-16-1.70002901e-16j,  2.91433544e-16+5.55111512e-17j,\n",
       "       -3.29597460e-17-8.32667268e-17j,  0.00000000e+00+6.93889390e-17j,\n",
       "        7.28583860e-17-2.22044605e-16j, -6.59194921e-17+5.55111512e-17j])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ifft(v) - inverse_fourier_transform(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Connection with the Laplacian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The functions $u_n$ (the columns of the Fourier transform matrix) are eigenvectors of the Laplacian:\n",
    "\n",
    "\\begin{align*}\n",
    "u_n(x)&:= \\cos\\left(2 \\pi \\frac{n}{N} x\\right) - i \\sin\\left(2 \\pi \\frac{n}{N} x\\right)\\\\\n",
    "\\Delta u_n(x)&:= \\left(-4 \\pi ^ 2\\frac{n^2}{N^2}\\right) u_n(x)\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "v * w \n",
    "= U^H ((U  w) \\odot (U  v))\n",
    "\\end{align*}\n",
    "\n",
    "or if $g_w=\\mbox{diag}(U w)$ is  filter\n",
    "\\begin{align*}\n",
    "v * w \n",
    "= U^H g_w U  w\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4.75972857, 3.84267706, 4.75547361, 4.69473343, 4.41980675,\n",
       "       4.38291796, 3.98819354, 4.36896188, 4.0910097 , 4.60256767,\n",
       "       4.9724941 , 5.12207561, 5.60658149, 4.99165311, 4.75559415,\n",
       "       4.99484192, 4.8482739 , 4.82087932, 4.5831123 , 4.63790021])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U = matrix_U(N)\n",
    "np.abs((U.conj().transpose() / N) @ ((U @ v) * (U @ w)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4.75972857, 3.84267706, 4.75547361, 4.69473343, 4.41980675,\n",
       "       4.38291796, 3.98819354, 4.36896188, 4.0910097 , 4.60256767,\n",
       "       4.9724941 , 5.12207561, 5.60658149, 4.99165311, 4.75559415,\n",
       "       4.99484192, 4.8482739 , 4.82087932, 4.5831123 , 4.63790021])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "conv(v, w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convolution on graphs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Plan**:\n",
    "    - Define the graph Laplacian\n",
    "    - Compute the spectrum\n",
    "    - Define a Fourier transform\n",
    "    - Define convolution on a graph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** From now on $G = (V, E)$ is an undirected, unweighted, simple graph."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Graph Laplacian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adjacency matrix\n",
    "\\begin{align*}\n",
    "A_{ij} = \\left\\{\n",
    "    \\begin{array}{ll}\n",
    "    1 &\\text{ if } e_{ij}\\in E\\\\\n",
    "    0 &\\text{ if } e_{ij}\\notin E\n",
    "    \\end{array}\n",
    "    \\right.\n",
    "\\end{align*}\n",
    "\n",
    "Degree matrix\n",
    "\\begin{align*}\n",
    "D_{ij} = \\left\\{\n",
    "    \\begin{array}{ll}\n",
    "    \\mbox{deg}(v_i) &\\text{ if } i=j\\\\\n",
    "    0 &\\text{ if } i\\neq j\n",
    "    \\end{array}\n",
    "    \\right.\n",
    "\\end{align*}\n",
    "\n",
    "Laplacian\n",
    "\\begin{align*}\n",
    "L &= D - A.\n",
    "\\end{align*}\n",
    "\n",
    "Normalized Laplacian\n",
    "\\begin{align*}\n",
    "L &= I - D^{-1/2} A D^{-1/2}.\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Graph spectrum, Fourier transform, and convolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Spectral decomposition of the Laplacian:\n",
    "\\begin{align*}\n",
    "L = U \\Lambda U^T\\\\\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "2. Fourier transform: if $v$ is a vector of features on the graph, then\n",
    "\\begin{align*}\n",
    "\\mathcal F (v) = U \\cdot v, \\;\\;\\mathcal F^{-1} (v) = U^T \\cdot v\\\\\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "3. Convolution with a filter $U \\cdot w$\n",
    "\\begin{align*}\n",
    "v * w = U ((U^T  w) \\odot (U^T  v) )\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "Or $g_w = \\mbox{diag}(U^T w)$ is a filter, then\n",
    "\\begin{align*}\n",
    "v * w = U g_w U^T  v\n",
    "\\end{align*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spectral-convolutional layers in PyTorch Geometric"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Problem:** Computing the spectrum is a global and very expensive property.\n",
    "\n",
    "**Goal:** Implementation as message passing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ChebConv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Original [paper](https://arxiv.org/pdf/1606.09375.pdf)\n",
    "- PyTorch [doc](https://pytorch-geometric.readthedocs.io/en/latest/modules/nn.html#torch_geometric.nn.conv.ChebConv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Goal: \n",
    "Compute $U g_w U^T x$ with $g_w = g_w(\\Lambda)$ a filter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Chebyshev approximation\n",
    "\n",
    "Chebyshev polynomials $T_k$:\n",
    "\\begin{align*}\n",
    "T_{k}(x) = 2 x T_{k-1}(x) - T_{k-2}(x), \\;\\; T_0(x) = 1, T_1(x) = x\n",
    "\\end{align*}\n",
    "\n",
    "#### Chebyshev approximation of the filter\n",
    "Aproximation of the filter:\n",
    "\\begin{align*}\n",
    "g_w(\\Lambda) = \\sum_{k=0}^K \\theta_k T_k(\\tilde \\Lambda),\\;\\;\\;\\;\\tilde \\Lambda = \\frac{2}{\\lambda_\\max} \\Lambda - I\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "#### Property\n",
    "If $L = U \\Lambda U^T$ then $T_k(L) = U T_k(\\Lambda) U^T$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fast approximated convolution \n",
    "\\begin{align*}\n",
    "v * w &= U g_w U^T x\n",
    "= U \\left(\\sum_{k=0}^K \\theta_k T_k(\\tilde \\Lambda) \\right)U^T x\n",
    "=\\sum_{k=0}^K  \\theta_k U  T_k(\\tilde \\Lambda) U^T x\\\\ \n",
    "&=\\sum_{k=0}^K  \\theta_k T_k(\\tilde L) x \n",
    "\\end{align*}\n",
    "\n",
    "\\begin{align*}\n",
    "\\tilde L = \\frac{2}{\\lambda_\\max} L - I\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Properties:\n",
    "- Depends on $L$ and $\\lambda_\\max$, not on $U, \\Sigma$\n",
    "- Uses only $K$-powers $\\Rightarrow$ only the $K$-th neighborhood of each node, localized filter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**As message passing:**\n",
    "![](fig/cheb_init.png)\n",
    "\n",
    "![](fig/cheb_norm.png)\n",
    "\n",
    "![](fig/cheb_forward.png)\n",
    "\n",
    "![](fig/cheb_message.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GCNConv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Original [paper](https://arxiv.org/pdf/1609.02907.pdf)\n",
    "- PyTorch [doc](https://pytorch-geometric.readthedocs.io/en/latest/modules/nn.html#torch_geometric.nn.conv.GCNConv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start from `ChebConv` and assume \n",
    "1. $K=1$ (linear approximation) so\n",
    "\\begin{align*}\n",
    "v * w \n",
    "&=\\sum_{k=0}^1  \\theta_k T_k(\\tilde L) x\n",
    "= \\theta_0 x + \\theta_1 \\tilde L x\\\\\n",
    "\\end{align*}\n",
    "\n",
    "2. $\\lambda_\\max =2$ so\n",
    "\\begin{align*}\n",
    "v * w \n",
    "&= \\theta_0 x + \\theta_1 (L - I) x\\\\\n",
    "&= \\theta_0 x - \\theta_1 D^{-1/2} A D^{1/2} x\\\\\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "3. $\\theta_0=-\\theta_1= \\theta$ so \n",
    "\\begin{align*}\n",
    "v * w = \\left(I + D^{-1/2} A D^{1/2}\\right) x \\theta\n",
    "\\end{align*}\n",
    "\n",
    "4. Renormalization of $\\theta$ by using \n",
    "\\begin{align*}\n",
    "\\tilde A&:= I + A\\\\\n",
    "\\tilde D_{ii}&:= \\sum_j \\tilde A_{ij}\n",
    "\\end{align*}\n",
    "so \n",
    "\\begin{align*}\n",
    "v * w = \\left(D^{-1/2} A D^{1/2}\\right) x \\theta\n",
    "\\end{align*}\n",
    "\n",
    "If $x$ is a $F$-dimensional feature vector, and we want an $F'$-dimensional feature vector as output:\n",
    "use $W'\\in \\mathbb{R}^{F\\times F'}$\n",
    "\\begin{align*}\n",
    "v * w = \\left(D^{-1/2} A D^{1/2}\\right) x \\Theta\n",
    "\\end{align*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nodewise:\n",
    "    ![image.png](fig/gcn_nodewise.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### As message passing\n",
    "See Tutorial3"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
